Comparing Age and Growth Increments from Bayesian and integrative data approaches for the deepwater snapper Pristipomoides filamentosus in the Hawaiian Islands

This Notebook contains code to run analysis and produce plots associated with the manuscript “Comparing Age and Growth Increments from Bayesian and integrative data approaches for the deepwater snapper Pristipomoides filamentosus in the Hawaiian Islands”

Workspace Setup

Clearing R enviornment

## Clearing R enviornment
rm(list = ls())
gc()

          used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
Ncells 2427637 129.7    4384820 234.2         NA  4384820 234.2
Vcells 4433993  33.9   10146329  77.5      16384  7317015  55.9

## Establishing Directory Heirarchy
proj_dir = getwd() 
# proj_dir = "/Volumes/GoogleDrive/My Drive/Weng Lab/Personal_Folders/Steve/dissertation work/Ch 4. Opakapaka Growth/Analysis"
data_dir = file.path(proj_dir, 'data')
src_dir = file.path(proj_dir, 'src')
results_dir = file.path(proj_dir, 'results')
run_results_dir = results_dir

set.seed(42)

Importing Packages and Dependencies

Loading saved workspace

After this script concludes, an workspace image will be saved so results can be accessed without re-running growth models. If this workspace already exists, we’ll load it in now.

if(file.exists(file.path(results_dir, 'Analysis Workspace.RData'))){
  # load(file.path(results_dir, 'Analysis Workspace.RData'))
}

User Defined Variables

min_time_at_lib = 60 # days

Bayesian Analysis (Models I - IV)

The following section uses JAGS to fit bayesian models under different constraints using the OTP data only.

Importing and preprocessing OTP Data

otp_data = load_okomoto_data()

NAs introduced by coercionNAs introduced by coercionNAs introduced by coercion

Formatting Data for Bayesian Modeling

The next sections create an list object named “data” consisting of the following objects:

L: A matrix of fork lengths (cm). Rows are individuals, columns are capture events dt: A matrix of delta t (years) corrosponding to time between capture events. n: A vector where values represent the number of valid recaptures for individuals N: A numeric variable corrosponding to the number of unique fish in our dataset

## Print number of rows in L. Should be 387
dim(L)[1] == 387

[1] TRUE

Specifying Bayesian Models

In this chunk, we define 4 models that differ in the way they that they treat VBGF parameters Linf and K.

Model 1: L infinity (Linf) and k will be allowed to vary between individuals. No Fixed Effects Model 2: Linf is allowed to vary between individuals. K is fixed Model 3: k is allowed to vary. Linf is fixed Model 4: Linf and k are fixed.

Note that some parameters in models 2 and 4 have been truncated to avoid upsetting the MCMC slicing algorthm

## Defining model 1
cat(
  '# Model 1
model{                                           
    for (i in 1:N)   {
        for (j in 2:n[i])   {
            L[i, j] ~ dnorm(L_Exp[i, j], tau)   
            L_Exp[i, j] <-  Linf[i] *(1.0 - exp(-k[i]*(A[i]+dt[i, j -1])))
            # posterior prediction
            L.pred[i, j] ~ dnorm(L_Exp[i, j], tau)
            p.value[i, j] <- step(L.pred[i, j] - L[i, j])
        }
        L[i, 1] ~ dnorm(L_Exp[i, 1], tau)
        L_Exp[i, 1] <-   Linf[i] *(1.0 - exp(-k[i]*A[i]))   

        # posterior prediction
        L.pred[i, 1] ~ dnorm(L_Exp[i, 1], tau)
        p.value[i, 1] <- step(L.pred[i, 1]- L[i, 1])

        Linf[i] ~ dnorm(Linf_mu,  Linf_tau)     
        k[i] ~ dnorm(k_mu, k_tau) T(0,1)
        A[i] ~ dgamma(Shape, rate)
    }
    Linf_std <- sqrt(1/Linf_tau)
    k_std <- sqrt(1/k_tau)
    variance <- 1/tau
    Linf_mu ~ dnorm(100, 0.0001)
    Linf_tau ~ dgamma(0.001, 0.0001)
    Shape ~ dunif(0, 100)
    rate ~ dunif(0, 100)
    k_mu ~ dbeta(1, 1)
    k_tau ~ dgamma(0.001, 0.0001)
    tau ~ dgamma(0.001, 0.0001)
}', 
file = file.path(src_dir, "VBGF JAGS Model 1.txt")
)

## Defining model 2
cat(
  '# Model 2
model{                                           
    for (i in 1:N)   {
        for (j in 2:n[i])   {
            L[i, j] ~ dnorm(L_Exp[i, j], tau)   
            L_Exp[i, j] <-  Linf[i] *(1.0 - exp(-k*(A[i]+dt[i, j -1])))
            L.pred[i, j] ~ dnorm(L_Exp[i, j], tau)
            p.value[i, j] <- step(L.pred[i, j] - L[i, j])
        }
        L[i, 1] ~ dnorm(L_Exp[i, 1], tau)
        L_Exp[i, 1] <-   Linf[i] *(1.0 - exp(-k*A[i]))  
        L.pred[i, 1] ~ dnorm(L_Exp[i, 1], tau)
        p.value[i, 1] <- step(L.pred[i, 1]- L[i, 1])
        Linf[i] ~ dnorm(Linf_mu,  Linf_tau)     
        A[i] ~ dgamma(Shape, rate)
    }
Linf_std <- sqrt(1/Linf_tau)
    k_std <- sqrt(1/k_tau)
    variance <- 1/tau
    k ~ dnorm(k_mu, k_tau) 
    Linf_mu ~ dnorm(100, 0.0001)
    Linf_tau ~ dgamma(0.001, 0.0001)
    Shape ~ dunif(0, 100)
    rate ~ dunif(0, 100)
    k_mu ~ dbeta(1, 1)  T(0.01,0.9)
    k_tau ~ dgamma(0.001 + 0.01, 0.0001) 
    tau ~ dgamma(0.001, 0.0001)

}
    ', 
file = file.path(src_dir, "VBGF JAGS Model 2.txt")
)

## Defining model 3
cat('
  # Model 3
model{                                           
    for (i in 1:N)   {
        for (j in 2:n[i])   {
            L[i, j] ~ dnorm(L_Exp[i, j], tau)   
            L_Exp[i, j] <-  Linf*(1.0 - exp(-k[i]*(A[i]+dt[i, j -1])))
            L.pred[i, j] ~ dnorm(L_Exp[i, j], tau)
            p.value[i, j] <- step(L.pred[i, j] - L[i, j])
        }
        L[i, 1] ~ dnorm(L_Exp[i, 1], tau)
        L_Exp[i, 1] <-   Linf *(1.0 - exp(-k[i]*A[i]))  
        L.pred[i, 1] ~ dnorm(L_Exp[i, 1], tau)
        p.value[i, 1] <- step(L.pred[i, 1]- L[i, 1])
        k[i] ~ dnorm(k_mu, k_tau) T(0,1)
        A[i] ~ dgamma(Shape, rate)
    }
    Linf_std <- sqrt(1/Linf_tau)
    k_std <- sqrt(1/k_tau)
    variance <- 1/tau
    Linf ~ dnorm(Linf_mu,  Linf_tau)
    Linf_mu ~ dnorm(100, 0.0001)
    Linf_tau ~ dgamma(0.01, 0.0001)
    Shape ~ dunif(0, 100)
    rate ~ dunif(0, 1000)
    k_mu ~ dbeta(1, 1)
    k_tau ~ dgamma(0.01, 0.0001)
    tau ~ dgamma(0.01, 0.0001)
}', 
file = file.path(src_dir, "VBGF JAGS Model 3.txt")
)

## Defining model 4
cat(
  '# Model 4
model{                                           
    for (i in 1:N)   {
        for (j in 2:n[i])   {
            L[i, j] ~ dnorm(L_Exp[i, j], tau)   
            L_Exp[i, j] <-  Linf*(1.0 - exp(-k*(A[i]+dt[i, j-1])))
            L.pred[i, j] ~ dnorm(L_Exp[i, j], tau)
            p.value[i, j] <- step(L.pred[i, j] - L[i, j])
        }
    # Predicting length at capture
        L[i, 1] ~ dnorm(L_Exp[i, 1], tau)
        L_Exp[i, 1] <-   Linf *(1.0 - exp(-k*A[i])) 
        L.pred[i, 1] ~ dnorm(L_Exp[i, 1], tau)
        p.value[i, 1] <- step(L.pred[i, 1]- L[i, 1])
        A[i] ~ dgamma(Shape, rate)
    }
    k_std <- sqrt(1/k_tau)
    variance <- 1/tau
    k ~ dnorm(k_mu, k_tau) 
    Linf ~ dnorm(Linf_mu,  Linf_tau)
    Linf_mu ~ dnorm(100, 0.0001)
    Linf_tau ~ dgamma(0.001, 0.0001)
    Linf_std <- sqrt(1/Linf_tau)
    Shape ~ dunif(0, 100)
    rate ~ dunif(0, 100)
    k_mu ~ dbeta(1, 1)  T(0.01,0.9)
    k_tau ~ dgamma(0.001, 0.0001) 
    tau ~ dgamma(0.001, 0.0001)
}
', 
file = file.path(src_dir, "VBGF JAGS Model 4.txt")
)

Initializing MCMC Chains

Here we’ll define inits, a list of initial starting points for our optimizer. inits is a list containing a set of lists corrosponding to each chain we’ll run in our optimizer. For this analyis there will be 3 chains. The first we’ll set using some reasonable estimates from our maximum likelihood work, the remaining 2 will use values 2X larger and smaller than the first chain.

## Initial values for each chain are stored in 3 lists with elements corrosponding to all variables to be initialized

Warning message:
In readChar(file, size, TRUE) : truncating string with embedded nuls

inits1 = list('Linf_mu' = 60,   'k_mu' = 0.30)
inits2 = list('Linf_mu' = 60*2, 'k_mu' = 0.3*2)
inits3 = list('Linf_mu' = 60/2, 'k_mu' = 0.3/2)
## Creating a single list that contains the lists corrosponding to each chain's initial values
inits = list(inits1, inits2, inits3)

## Chain parameters
# The number of samples in the posterior distribution = (n_iterations - n_burnin) / n_thin 
n_iterations = 200# 500000 # How many total iterations to run
n_burnin = 20 # 10000 # How many iterations to run during the model's burn in phase 
n_thin = 2 # 50 # Retain one in n_thin samples from the posterior distribution. 

Running Bayesian Models

## Model 1: Linf and K vary between individuals
model_1 = jags(data, inits, 
               model.file = file.path(src_dir, "VBGF JAGS Model 1.txt"),
               parameters.to.save =  c('Linf_mu', 'Linf_std', 'Linf_tau', 'Shape', 'k_mu', 'k_std', 'k_tau', 'rate', 'tau', 'variance'),
               DIC = T, 
               n.chains = 3, n.iter = n_iterations, n.burnin = n_burnin, n.thin = n_thin)

module glm loaded

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 803
   Unobserved stochastic nodes: 1971
   Total graph size: 9210

Initializing model


  |                                                        
  |                                                  |   0%
  |                                                        
  |*                                                 |   2%
  |                                                        
  |**                                                |   4%
  |                                                        
  |***                                               |   7%
  |                                                        
  |****                                              |   9%
  |                                                        
  |******                                            |  11%
  |                                                        
  |*******                                           |  13%
  |                                                        
  |********                                          |  16%
  |                                                        
  |*********                                         |  18%
  |                                                        
  |**********                                        |  20%
  |                                                        
  |***********                                       |  22%
  |                                                        
  |************                                      |  24%
  |                                                        
  |*************                                     |  27%
  |                                                        
  |**************                                    |  29%
  |                                                        
  |****************                                  |  31%
  |                                                        
  |*****************                                 |  33%
  |                                                        
  |******************                                |  36%
  |                                                        
  |*******************                               |  38%
  |                                                        
  |********************                              |  40%
  |                                                        
  |*********************                             |  42%
  |                                                        
  |**********************                            |  44%
  |                                                        
  |***********************                           |  47%
  |                                                        
  |************************                          |  49%
  |                                                        
  |**************************                        |  51%
  |                                                        
  |***************************                       |  53%
  |                                                        
  |****************************                      |  56%
  |                                                        
  |*****************************                     |  58%
  |                                                        
  |******************************                    |  60%
  |                                                        
  |*******************************                   |  62%
  |                                                        
  |********************************                  |  64%
  |                                                        
  |*********************************                 |  67%
  |                                                        
  |**********************************                |  69%
  |                                                        
  |************************************              |  71%
  |                                                        
  |*************************************             |  73%
  |                                                        
  |**************************************            |  76%
  |                                                        
  |***************************************           |  78%
  |                                                        
  |****************************************          |  80%
  |                                                        
  |*****************************************         |  82%
  |                                                        
  |******************************************        |  84%
  |                                                        
  |*******************************************       |  87%
  |                                                        
  |********************************************      |  89%
  |                                                        
  |**********************************************    |  91%
  |                                                        
  |***********************************************   |  93%
  |                                                        
  |************************************************  |  96%
  |                                                        
  |************************************************* |  98%
  |                                                        
  |**************************************************| 100%

## Model 2: Linf varies between individuals, K is fixed
model_2 = jags(data, inits, 
               model.file = file.path(src_dir, "VBGF JAGS Model 2.txt"),
               parameters.to.save =  c('Linf_mu', 'Linf_std', 'Linf_tau', 'Shape', 'k_mu', 'k_std', 'k_tau', 'rate', 'tau', 'variance'),
               DIC = T, 
               n.chains = 3, n.iter = n_iterations, n.burnin = n_burnin, n.thin = n_thin)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 803
   Unobserved stochastic nodes: 1585
   Total graph size: 8441

Initializing model


  |                                                        
  |                                                  |   0%
  |                                                        
  |*                                                 |   2%
  |                                                        
  |**                                                |   4%
  |                                                        
  |***                                               |   7%
  |                                                        
  |****                                              |   9%
  |                                                        
  |******                                            |  11%
  |                                                        
  |*******                                           |  13%
  |                                                        
  |********                                          |  16%
  |                                                        
  |*********                                         |  18%
  |                                                        
  |**********                                        |  20%
  |                                                        
  |***********                                       |  22%
  |                                                        
  |************                                      |  24%
  |                                                        
  |*************                                     |  27%
  |                                                        
  |**************                                    |  29%
  |                                                        
  |****************                                  |  31%
  |                                                        
  |*****************                                 |  33%
  |                                                        
  |******************                                |  36%
  |                                                        
  |*******************                               |  38%
  |                                                        
  |********************                              |  40%
  |                                                        
  |*********************                             |  42%
  |                                                        
  |**********************                            |  44%
  |                                                        
  |***********************                           |  47%
  |                                                        
  |************************                          |  49%
  |                                                        
  |**************************                        |  51%
  |                                                        
  |***************************                       |  53%
  |                                                        
  |****************************                      |  56%
  |                                                        
  |*****************************                     |  58%
  |                                                        
  |******************************                    |  60%
  |                                                        
  |*******************************                   |  62%
  |                                                        
  |********************************                  |  64%
  |                                                        
  |*********************************                 |  67%
  |                                                        
  |**********************************                |  69%
  |                                                        
  |************************************              |  71%
  |                                                        
  |*************************************             |  73%
  |                                                        
  |**************************************            |  76%
  |                                                        
  |***************************************           |  78%
  |                                                        
  |****************************************          |  80%
  |                                                        
  |*****************************************         |  82%
  |                                                        
  |******************************************        |  84%
  |                                                        
  |*******************************************       |  87%
  |                                                        
  |********************************************      |  89%
  |                                                        
  |**********************************************    |  91%
  |                                                        
  |***********************************************   |  93%
  |                                                        
  |************************************************  |  96%
  |                                                        
  |************************************************* |  98%
  |                                                        
  |**************************************************| 100%

## Model 3: K varies between individuals, Linf is fixed
model_3 = jags(data, inits, 
               model.file = file.path(src_dir, "VBGF JAGS Model 3.txt"),
               parameters.to.save =  c('Linf_mu', 'Linf_std', 'Linf_tau', 'Shape', 'k_mu', 'k_std', 'k_tau', 'rate', 'tau', 'variance'),
               DIC = T, 
               n.chains = 3, n.iter = n_iterations, n.burnin = n_burnin, n.thin = n_thin)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 803
   Unobserved stochastic nodes: 1585
   Total graph size: 8825

Initializing model


  |                                                        
  |                                                  |   0%
  |                                                        
  |*                                                 |   2%
  |                                                        
  |**                                                |   4%
  |                                                        
  |***                                               |   7%
  |                                                        
  |****                                              |   9%
  |                                                        
  |******                                            |  11%
  |                                                        
  |*******                                           |  13%
  |                                                        
  |********                                          |  16%
  |                                                        
  |*********                                         |  18%
  |                                                        
  |**********                                        |  20%
  |                                                        
  |***********                                       |  22%
  |                                                        
  |************                                      |  24%
  |                                                        
  |*************                                     |  27%
  |                                                        
  |**************                                    |  29%
  |                                                        
  |****************                                  |  31%
  |                                                        
  |*****************                                 |  33%
  |                                                        
  |******************                                |  36%
  |                                                        
  |*******************                               |  38%
  |                                                        
  |********************                              |  40%
  |                                                        
  |*********************                             |  42%
  |                                                        
  |**********************                            |  44%
  |                                                        
  |***********************                           |  47%
  |                                                        
  |************************                          |  49%
  |                                                        
  |**************************                        |  51%
  |                                                        
  |***************************                       |  53%
  |                                                        
  |****************************                      |  56%
  |                                                        
  |*****************************                     |  58%
  |                                                        
  |******************************                    |  60%
  |                                                        
  |*******************************                   |  62%
  |                                                        
  |********************************                  |  64%
  |                                                        
  |*********************************                 |  67%
  |                                                        
  |**********************************                |  69%
  |                                                        
  |************************************              |  71%
  |                                                        
  |*************************************             |  73%
  |                                                        
  |**************************************            |  76%
  |                                                        
  |***************************************           |  78%
  |                                                        
  |****************************************          |  80%
  |                                                        
  |*****************************************         |  82%
  |                                                        
  |******************************************        |  84%
  |                                                        
  |*******************************************       |  87%
  |                                                        
  |********************************************      |  89%
  |                                                        
  |**********************************************    |  91%
  |                                                        
  |***********************************************   |  93%
  |                                                        
  |************************************************  |  96%
  |                                                        
  |************************************************* |  98%
  |                                                        
  |**************************************************| 100%

## Model 4: Linf and K are fixed
model_4 = jags(data, inits, 
               model.file = file.path(src_dir, "VBGF JAGS Model 4.txt"),
               parameters.to.save =  c('Linf_mu', 'Linf_std', 'Linf_tau', 'Shape', 'k_mu', 'k_std', 'k_tau', 'rate', 'tau', 'variance'),
               DIC = T, 
               n.chains = 3, n.iter = n_iterations, n.burnin = n_burnin, n.thin = n_thin)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 803
   Unobserved stochastic nodes: 1199
   Total graph size: 8054

Initializing model


  |                                                        
  |                                                  |   0%
  |                                                        
  |*                                                 |   2%
  |                                                        
  |**                                                |   4%
  |                                                        
  |***                                               |   7%
  |                                                        
  |****                                              |   9%
  |                                                        
  |******                                            |  11%
  |                                                        
  |*******                                           |  13%
  |                                                        
  |********                                          |  16%
  |                                                        
  |*********                                         |  18%
  |                                                        
  |**********                                        |  20%
  |                                                        
  |***********                                       |  22%
  |                                                        
  |************                                      |  24%
  |                                                        
  |*************                                     |  27%
  |                                                        
  |**************                                    |  29%
  |                                                        
  |****************                                  |  31%
  |                                                        
  |*****************                                 |  33%
  |                                                        
  |******************                                |  36%
  |                                                        
  |*******************                               |  38%
  |                                                        
  |********************                              |  40%
  |                                                        
  |*********************                             |  42%
  |                                                        
  |**********************                            |  44%
  |                                                        
  |***********************                           |  47%
  |                                                        
  |************************                          |  49%
  |                                                        
  |**************************                        |  51%
  |                                                        
  |***************************                       |  53%
  |                                                        
  |****************************                      |  56%
  |                                                        
  |*****************************                     |  58%
  |                                                        
  |******************************                    |  60%
  |                                                        
  |*******************************                   |  62%
  |                                                        
  |********************************                  |  64%
  |                                                        
  |*********************************                 |  67%
  |                                                        
  |**********************************                |  69%
  |                                                        
  |************************************              |  71%
  |                                                        
  |*************************************             |  73%
  |                                                        
  |**************************************            |  76%
  |                                                        
  |***************************************           |  78%
  |                                                        
  |****************************************          |  80%
  |                                                        
  |*****************************************         |  82%
  |                                                        
  |******************************************        |  84%
  |                                                        
  |*******************************************       |  87%
  |                                                        
  |********************************************      |  89%
  |                                                        
  |**********************************************    |  91%
  |                                                        
  |***********************************************   |  93%
  |                                                        
  |************************************************  |  96%
  |                                                        
  |************************************************* |  98%
  |                                                        
  |**************************************************| 100%

growth_models = list('model_1' = model_1, 'model_2' = model_2, 'model_3' = model_3, 'model_4' = model_4)

Model Diagnostics

Producing diagnostic plots and checking convergence for each of our 4 growth models,

##### Model Diagnostics
for(i in 1:length(growth_models)){
  # readline(paste('Press Enter to View Diagnostics for', names(growth_models)[i]))
  model = growth_models[[i]]
  print(paste("Diagnostic plots for", names(growth_models)[i]))
  
  ## Write out model summary table
  # write.csv(model$BUGSoutput[10], file.path(results_dir, paste( names(growth_models)[i], ' Parameter Summaries.csv')))
}

[1] "Diagnostic plots for model_1"
[1] "Diagnostic plots for model_2"
[1] "Diagnostic plots for model_3"
[1] "Diagnostic plots for model_4"

for(i in 1:length(growth_models)){
  ### Summary Statistics and Parameter Estimates
  summary(model)
  
  ### Generating MCMC object for analysis
  model_mcmc = as.mcmc(model)
  
  par(mfrow = c(2, 4))
  ### xy plot
  xyplot(model_mcmc)
  
  ### Trace and Density Plots
  plot(model_mcmc)
  
  ### Autocorrelation plots
  autocorr.plot(model_mcmc)
  
  #### Other Diagnostics using CODA package
  gelman.plot(model_mcmc)
  geweke.diag(model_mcmc)
  geweke.plot(model_mcmc)
  raftery.diag(model_mcmc)
  heidel.diag(model_mcmc)
}

dev.off()

null device 
          1 

Our diagnostic plots don’t look bad and we can see that the Gelman Rubin convergene plots indicate that all models converged within the model burn-in phase.

After reviewing our diagnostic plots, we’ll move on to looking at which term in our models have the most credibility. We’ll do this using the coefficient of variation, equivilant to the mean divided by the standard deviation. A high coefficient of variability is indicitive of a poor parameter fit.

We assume the model 1, which allows both L and K to be fit independently for each fish is the best model. A low coefficient of variation for parameters under this model should confirm this. We’ll then look at the remaining 3 models to see how variability in model terms is affected by fixing a given parameter across the population.

First we need to calculate the coefficient of variation for Linf_mu and k_mu parameters for each model. Then we’ll make a plot examining the magnitude of the coefficient of variation for each model’s parameters.

#### Extracting coefficients of variation for Linf and K paramters
cv_df = data.frame(stringsAsFactors = F)
for(i in 1:length(growth_models)){
  curr_mod = growth_models[[i]]
  cv_df = rbind(cv_df, data.frame('model' = names(growth_models)[i], 'Parameter' = 'Linf', 'cv' = curr_mod$BUGSoutput$summary["Linf_mu",'sd'] / curr_mod$BUGSoutput$summary["Linf_mu","mean"] * 100), 
                data.frame('model' = names(growth_models)[i], 'Parameter' = 'k', 'cv' = curr_mod$BUGSoutput$summary["k_mu",'sd'] / curr_mod$BUGSoutput$summary["k_mu","mean"] * 100))
}

## Adding the source of variability for each model and term to our data frame
cv_df$`source of individual variability` = 'Other'
cv_df$`source of individual variability`[cv_df$model == 'model_1'] = 'Both'
cv_df$`source of individual variability`[cv_df$model == 'model_4'] = 'Neither'
cv_df$`source of individual variability`[cv_df$model == 'model_2' & cv_df$Parameter == 'Linf'] = 'Self'
cv_df$`source of individual variability`[cv_df$model == 'model_3' & cv_df$Parameter == 'k'] = 'Self'
cv_df$`source of individual variability` = as.factor(cv_df$`source of individual variability`)
cv_df$`source of individual variability` = fct_relevel(cv_df$`source of individual variability`, c('Both', 'Self', 'Other', 'Neither'))

Visualizing CV

### Plotting our coefficient of variation
fig3 = ggplot(
              data = cv_df, 
              mapping = aes(x = `source of individual variability`, y = cv, col = Parameter)
              ) + 
        geom_point() + 
        geom_line(aes(group = Parameter)) + 
        labs(
          x = 'Source of Individual Variability', 
          y = 'Coefficient of Variation (Percent)', 
          fill = 'Parameter'
          ) + 
        theme(
          legend.justification = c(1, 1), 
          legend.position = c(.15, .95)
          ) + 
        ylim(0,100)

## Save Figure
  # ggsave(
  #   filename = 'Fig3 - Coefficients of Variation.pdf', 
  #   plot = fig3,
  #   device = 'pdf', 
  #   path = results_dir
  #   )
  
## Show figure as output
  print(fig3)

The x-axis labels indicate wether a term is fixed or not. “Both” refers to model 1 as both terms are fit independently. “Self”" and “Other” are models 2 and 3, while “neither” is model 4. We see that when both parameters are fit independently (model 1), both terms have the lowest coefficient of variation, indicating it is the model with the best fit. We see, unsuprisingly, as we make our effects fixed across the population, the deviation associated increases.

While model 1 is clearly the best performing model, when we compare parameter estimates from Models 1 and 2, we can infer the model 2, is likely credible as well.

The additional Models 2-4 suggested that individual variability in both K and L_∞ was important, with perhaps variability in L_∞ being more important based upon the response of L_∞ standard deviation from the base case of Model 1 to the constrained individual variability in Model 3 and Model 4.

Maximum Likelihood Integrative Models

Formatting Tagging Data for ML integrative models

Now we want to format a table with lm (length at marking), lr (length at recapture), and dt (elapsed time) Note: If a fish was recaptured multiple times, there is a single entry for that individual corrosponding to the length at initial marking and the length at final recapture

tagdat = data.frame(
            'L1' = L[ ,1],
            "L2" = rep(0, length(L[,1])), 
            " "  = rep(0, length(L[,1])),  # A dummy column. This is required for matrix indexing in Laslett's utility functions
            "dt" = rep(0, length(L[,1])),
            "L2measurer" = rep(0, length(L[,1]))
            )

for (i in 1:nrow(tagdat)){
  tagdat$L2[i] = L[i, max(which(!is.na(L[i, ])))]
  tagdat$dt[i] = sum(dt[i], na.rm = T)
}

Importing Additional Data Sets - Otolith Data and Length Frequency Data

#### Otolith Data (Ralston and Miyamoto 1983, DeMartini et al. 1994, Andrews et al. 2012) 
otodat = read.csv(file.path(data_dir, "RalstonMiyamotoandDemartiniAndrews.csv"), col.names = c("age", "len", "source"))

length_frequency_data = read.csv(file.path(data_dir, 'moffit and parrish pseudo lenght frequency data.csv'))
  length_frequency_data$date = as.POSIXct(length_frequency_data$date)

Decomposing length Frequency data

Performing modal length frequency decomposition tracking maximum monthly mode fork lengths. (See Eveson 2014? for details) Fork lengths are related to age based on spawning time and time to recruitment

lfdat = length_freq_decomp(length_frequency_data, plot = TRUE, fixed_modes = TRUE)

number of iterations= 21 
number of iterations= 26 
number of iterations= 25 
number of iterations= 72 
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
number of iterations= 1 
number of iterations= 12 
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
One of the variances is going to zero;  trying new starting values.
number of iterations= 1 
number of iterations= 29 

Fitting Preliminary Models

Using the methods of Laslett et al, fit separate models to tagging data, length-age data, and length frequency data

#### Initializing global model parameters for all ML models
growth.ssnl.f<- growthvb.f
npf <- 1  #number of parameters passed to growth.ssnl.f (in this case k)
npA <- 2  #number of parameters in distribution of release ages for tag model

#Specifying starting parameters, as well as upper and lower bounds for parameter estimation
##       mu.L, sig.L,  k,  mu.A, sig.A, sig.sci, sig.f, a0, sig.oto, sig.lf
p0 <- c(  70,     1, .10,     1,   .10,      1,     0,   0,     1,     0)
lb <- c(  50,   0.1, .05,   0.1,   .05,    0.1,     0,   0,     0,     0)
ub <- c( 110,  15.0, .50,   1.5,   .50,   15.0,     0,   0,     0,     0)

### 1. Mark Recapture Data
vbgf_growth_increment <- nlminb(p0,joint.logl.f,lower=lb,upper=ub,npf=npf,npA=npA,tagdat=tagdat,otodat=otodat,lfdat=lfdat, wt.oto=0,wt.tag=1,wt.lf=0)

[1] "maximum on upper bound"
[1] "rows"
  [1]  22  23  24  26  27  28  29  30  32  34  35  36  37  39  40  41  42  43  44  50  51  52  54  55  56  58  59
 [28]  60  61  63  64  66  77  78  79  80  85  86  89  90  94  99 101 102 111 113 114 115 121 124 129 130 134 140
 [55] 141 143 144 145 148 152 154 157 158 159 163 178 179 182 184 185 189 191 193 197 203 207 209 211 212 214 218
 [82] 224 225 226 228 229 232 237 238 247 252 253 266 267 269 270 271 272 276 277 278 279 283 287 288 300 301 306
[109] 308 310 311 312 313 314 315 316 317 318 319 320 324 330 333 336 341 353 354 355 356 357 359 361 362 364 365
[136] 366 367 368 370 371 372 373 377 380 381 382 384

Error in lna[tf] <- log(a[tf]) : 
  NAs are not allowed in subscripted assignments

vbgf_length_age

$par
 [1] 63.6466006  6.9670575  0.3680429  1.0000000  0.1000000  1.0000000  0.0000000  0.0000000  0.8463622
[10]  1.0000000

$objective
[1] 381.5317

$convergence
[1] 0

$iterations
[1] 41

$evaluations
function gradient 
      48      430 

$message
[1] "relative convergence (4)"

vbgf_length_frequency_unconstrained

$par
 [1] 120.646746   1.000000   0.050000   1.000000   0.100000   1.000000   0.000000  -1.761237   1.000000
[10]   0.791511

$objective
[1] 26.6886

$convergence
[1] 0

$iterations
[1] 140

$evaluations
function gradient 
     187     1420 

$message
[1] "relative convergence (4)"

### Saving best growth parameter estimates from each model to a table
results = data.frame(stringsAsFactors = FALSE)
  results = rbind(results, cbind('Model 5 - Mark Recapture - All Data', t(as.vector(vbgf_growth_increment$par))))
  results = rbind(results, cbind('Length at Age', t(as.vector(vbgf_length_age$par))))
  results = rbind(results, cbind('Length Frequency (Unconstrained)', t(as.vector(vbgf_length_frequency_unconstrained$par))))
  results = rbind(results, cbind('Length Frequency (Constrained)', t(as.vector(vbgf_length_frequency_constrained$par))))

Integrative Model Fitting

### Setting intial params for Integrative Models
##       mu.L, sig.L,  k,  mu.A, sig.A, sig.sci, sig.f, a0, sig.oto, sig.lf
p0 <- c(  70,     1, .10,     1,   .10,      1,     0,   0,     1,      1)
lb <- c(  40,   0.01, .05,   0.1,   .05,    0.1,     0,  -10,  0.1,    0.1)
ub <- c( 110,  15.0, .50,   1.5,   .50,   15.0,     0,   10,   15,     15)

Fitting Integrative Models

### 6. Model including all Data sources - Equal weighting to each data type
fit.vb.equalwt.grouped <- nlminb(p0, joint.logl.f, lower=lb, upper=ub, npf=npf, npA=npA, tagdat=tagdat, otodat=otodat, lfdat=lfdat, wt.oto=1/dim(otodat)[1], wt.tag=1/dim(tagdat)[1], wt.lf=1/length(lfdat$curr_month_year))
results = rbind(results, cbind('Model 6 - All Data - Equal Weighting', t(as.vector(fit.vb.equalwt.grouped$par))))

### 7. Model including all Data sources - weighting based on number of sample size
fit.vb.byn.grouped <- nlminb(p0, joint.logl.f, lower=lb, upper=ub, npf=npf, npA=npA, tagdat=tagdat, otodat=otodat, lfdat=lfdat, wt.oto=1, wt.tag=1, wt.lf=1)
results = rbind(results, cbind('Model 7 - All Data - Weighted by n', t(as.vector(fit.vb.byn.grouped$par))))

### 8. Model including all Data sources treated individually - with equal weighting
fit.vb.equalwt.indv <- nlminb(p0, joint.logl.f, lower=lb, upper=ub, npf=npf, npA=npA, tagdat=tagdat, tagdat2 = NULL, otodat=otodat[otodat$source == 'demartini', ], otodat2=otodat[otodat$source == 'ralston and miyamoto', ], otodat3=otodat[otodat$source == 'andrews bomb carbon', ], otodat4=otodat[otodat$source == 'andrews lead radium', ], lfdat=lfdat, lfdat2=NULL, wt.oto= 1/dim(otodat[otodat$source == 'demartini', ])[1], wt.oto2= 1/dim(otodat[otodat$source == 'ralston and miyamoto', ])[1], wt.oto3=1/dim(otodat[otodat$source == 'andrews bomb carbon', ])[1], wt.oto4=1/dim(otodat[otodat$source == 'andrews lead radium', ])[1], wt.tag = 1/dim(tagdat)[1], wt.tag2 = 0, wt.lf = 1/length(lfdat$curr_month_year), wt.lf2 = 0)
results = rbind(results, cbind('Model 8 - Separated Data - Equal Weighting', t(as.vector(fit.vb.equalwt.indv$par))))

### 9. Model including all Data sources treated individually - weighting based on number of sample size
fit.vb.byn.indv <- nlminb(p0, joint.logl.f, lower=lb, upper=ub, npf=npf, npA=npA, tagdat=tagdat, tagdat2 = NULL, otodat=otodat[otodat$source == 'demartini', ], otodat2=otodat[otodat$source == 'ralston and miyamoto', ], otodat3=otodat[otodat$source == 'andrews bomb carbon', ], otodat4=otodat[otodat$source == 'andrews lead radium', ], lfdat=lfdat, lfdat2=NULL, wt.oto= 1, wt.oto2= 1, wt.oto3=1, wt.oto4=1, wt.tag = 1, wt.tag2 = 0, wt.lf = 1, wt.lf2 = 0)
results = rbind(results, cbind('Model 9 - Separated Data - Weighted by n', t(as.vector(fit.vb.byn.indv$par))))

### 10. Model without Ralston & Miyamoto - Equal weighting (Because Brett said this was shit!)
fit.vb.byn.indv.no.ralston <- nlminb(p0, joint.logl.f, lower=lb, upper=ub, npf=npf, npA=npA, tagdat=tagdat, tagdat2 = NULL, otodat=otodat[otodat$source == 'demartini', ], otodat2=otodat[otodat$source == 'ralston and miyamoto', ], otodat3=otodat[otodat$source == 'andrews bomb carbon', ], otodat4=otodat[otodat$source == 'andrews lead radium', ], lfdat=lfdat, lfdat2=NULL, wt.oto= 1/dim(otodat[otodat$source == 'demartini', ])[1], wt.oto2= 0, wt.oto3=1/dim(otodat[otodat$source == 'andrews bomb carbon', ])[1], wt.oto4=1/dim(otodat[otodat$source == 'andrews lead radium', ])[1], wt.tag = 1/dim(tagdat)[1], wt.tag2 = 0, wt.lf = 1/length(lfdat$curr_month_year), wt.lf2 = 0)
results = rbind(results, cbind('Model 10 - Separated Data - Equal Weighting - No R&M', t(as.vector(fit.vb.byn.indv.no.ralston$par))))

### 11. Model without Ralston & Miyamoto - weighted by n (Because Brett said this was shit!)
fit.vb.byn.indv.no.ralston <- nlminb(p0, joint.logl.f, lower=lb, upper=ub, npf=npf, npA=npA, tagdat=tagdat, tagdat2 = NULL, otodat=otodat[otodat$source == 'demartini', ], otodat2=otodat[otodat$source == 'ralston and miyamoto', ], otodat3=otodat[otodat$source == 'andrews bomb carbon', ], otodat4=otodat[otodat$source == 'andrews lead radium', ], lfdat=lfdat, lfdat2=NULL, wt.oto= 1, wt.oto2= 0, wt.oto3=1, wt.oto4=1, wt.tag = 1, wt.tag2 = 0, wt.lf = 1, wt.lf2 = 0)
results = rbind(results, cbind('Model 11 - Separated Data - Weighted by n - No R&M', t(as.vector(fit.vb.byn.indv.no.ralston$par))))
``

Error: attempt to use zero-length variable name

Comparing Model Structures to one another and to literature values

Determining preferred model structure using a train/test split bootstrapping procedure with 2/3 of the data allocated to model training and the remaining 1/3 of the data used to calculate RMSE.

These results will be compared to existing literature values and our estimates from Models 1-4. Here these values are imported as a .csv file

lit_vbgc_params = read.csv(file.path(data_dir, "Parameter Estimates.csv"), stringsAsFactors = FALSE)
lit_vbgc_params = lit_vbgc_params[!is.na(lit_vbgc_params$Linf), ]
colnames(lit_vbgc_params) = c('author', 'n', 'linf', 'k', 't0', 'region', 'method')
lit_vbgc_params = lit_vbgc_params[c(1:20, 22:25), ]

Iterative fitting and model evaluation using RMSE

n_train = round(dim(tagdat)[1] * (2/3))
n_iterations = 10000

mod_eval_results = evaluate_models(cross_validation_iterations = n_iterations)
print(mod_eval_results)

Visualizing model comparison

## Reshaping data to make boxplots
mod_eval_results_lf = as.data.frame(t(mod_eval_results[ ,1:7]))
mod_eval_results_lf$model_id = rownames(mod_eval_results_lf)
mod_eval_results_lf = reshape(mod_eval_results_lf, varying = colnames(mod_eval_results_lf[1:n_iterations]), idvar = 'model_id', direction = "long")
mod_eval_results_lf = mod_eval_results_lf

### Visualizing Results - Boxplot
boxplot(mod_eval_results_lf$result ~ mod_eval_results_lf$model_id, ylim = c(0, 7))

Declairing prefered integrative model structure

integrative_models = mod_eval_results[ ,2:7]
integrative_model_scores = c()
for(i in 1:dim(integrative_models)[1]){
  integrative_model_scores = c(integrative_model_scores, names(which.min(integrative_models[i, ])))
}
## Which model was most frequently the best one?
best_integrative_model = names(which.max(table(integrative_model_scores)))

Cleaning up results summary and writing it out to a .csv

colnames(results) = c('model_id', 'mu.L', 'sig.L',  'k',  'mu.A', 'sig.A', 'sig.sci', 'sig.f', 'a0', 'sig.oto', 'sig.lf')
results$`time to 90%` = yrs_to_.9_linf(linf = as.numeric(results$mu.L), k = as.numeric(results$k), a0 = as.numeric(results$a0))
# write.csv(results, file = file.path(run_results_dir, 'likelihood_parameter_estimates_with_full_data.csv'))

Comparing Best Model to Model without additional datasources (Model 5)

## Now comparing best integrative model to just tagging data
integrative_vs_tagging = mod_eval_results[ ,which(colnames(mod_eval_results) %in% c('model 5', best_integrative_model))]
integrative_vs_tagging_model_scores = c()
for(i in 1:dim(integrative_models)[1]){
  integrative_vs_tagging_model_scores = c(integrative_vs_tagging_model_scores, names(which.min(integrative_vs_tagging[i, ])))
}
best_model = names(which.max(table(integrative_vs_tagging_model_scores)))

### visualizing results
pdf(file.path(run_results_dir, 'Barplot of tagging vs. best integrative model.pdf'), width = 11, height = 8.5)
  barplot(prop.table(table(integrative_vs_tagging_model_scores)))
dev.off()

Comparing best model against all prviously published parameters for the region

nll_names = colnames(mod_eval_results)[1:7]
lit_names = colnames(mod_eval_results)[8:18]
bayes_names = colnames(mod_eval_results)[19:22]
lit_vs_int_vs_bayes = mod_eval_results[ ,colnames(mod_eval_results) %in% c(best_integrative_model, lit_names)]
lit_vs_int_vs_bayes_scores = c()
for(i in 1:dim(integrative_models)[1]){
  lit_vs_int_vs_bayes_scores = c(lit_vs_int_vs_bayes_scores, names(which.min(lit_vs_int_vs_bayes[i, ])))
}
best_model_lit_bayes_integrated = names(which.max(table(lit_vs_int_vs_bayes_scores)))

pdf(file.path(run_results_dir, 'Barplot of lit vs. bayes vs. best integrative model.pdf'), width = 11, height = 8.5)
par(las = 2)  
barplot(prop.table(table(lit_vs_int_vs_bayes_scores)))
dev.off()

### Subsetting model structures 6-11
nll_eval_results = mod_eval_results[, 2:7]

## Determining the number of NA iterations
na_index = c()
for(i in 1:length(nll_eval_results[ ,1])){
  if(any(is.na(nll_eval_results[i, ]))){
    na_index = c(na_index, i)
  }
}
na_index = unique(na_index)
### How many iterations failed to converge?
print(paste('Iterations failing to converge:', length(na_index)))
# nll_eval_results = nll_eval_results[-na_index, ]

## Getting summary stats for NLL models
print('Summary stats of competing model structures')
print(paste('Range: ', round(range(nll_eval_results, na.rm = TRUE), 2)))
nll_vec = as.vector(nll_eval_results)
nll_vec = nll_vec[!is.na(nll_vec)]
print(paste('mean:', round(mean(nll_vec), 2)))
print(paste('standard deviation:', round(sd(nll_vec), 2)))

### Determining which model performed best over all iterations
best_models = c()
for(i in 1:dim(nll_eval_results)[1]){
  best_models = c(best_models, names(which.min(nll_eval_results[i, ])))
}
print('Best Models')
table(best_models)

### Getting stats on the best performing model
print('Summary Stats for prefered integrated model')
best_nll_mod = mod_eval_results[ ,best_model]
print(paste('range:', round(range(best_nll_mod, na.rm = TRUE), 2)))
print(paste('mean:', round(mean(best_nll_mod, na.rm = TRUE), 2)))
print(paste('standard deviation:', round(sd(best_nll_mod, na.rm = TRUE), 2)))

### Getting stats on the model based only on tagging data
print('Summary Stats for Tagging Only Model (Model 5)')
mod_5 = as.vector(mod_eval_results[ ,'model 5'])
print(paste('range:', round(range(mod_5, na.rm = TRUE), 2)))
print(paste('mean:', round(mean(mod_5[!is.na(mod_5)]), 2)))
print(paste('standard deviation:', round(sd(mod_5[!is.na(mod_5)]), 2)))

### Comparing the perfered model to the tagging data only model
print('Comparing prefered integrative and tagging only models')
tagging_vs_integrative_df = cbind(mod_eval_results$`model 5`, mod_eval_results[ ,best_model])
colnames(tagging_vs_integrative_df) = c('model 5', best_model)

tagging_vs_integrative = c()
for(i in 1:length(tagging_vs_integrative_df[ ,1])){
  tagging_vs_integrative = c(tagging_vs_integrative, colnames(tagging_vs_integrative_df)[which.min(tagging_vs_integrative_df[i, ])])
}
table(tagging_vs_integrative)

### Summary stats on tagging and integrative models
pred_var_diff_tvc = tagging_vs_composite_df[ ,1] - tagging_vs_composite_df[ ,2]
print(paste('range in predicteve difference:', round(range(pred_var_diff_tvc, na.rm = TRUE), )))
print(paste('mean:', round(mean(pred_var_diff_tvc, na.rm = TRUE), 2)))
print(paste('standard deviation:', round(sd( as.vector(pred_var_diff_tvc)[!is.na(as.vector(pred_var_diff_tvc))]), 2)))

#### Getting summary stats on all literature models 
print('Summary Statistics for Literature Models')
lit = mod_eval_results[, 8:18]
print(paste('range:', round(range(lit, na.rm = TRUE),2)))
lit_vec = as.vector(lit)
lit_vec = lit_vec[!is.na(lit_vec)]
print(paste('mean:', round(mean(lit_vec), 2)))
print(paste('Standard Deviation:', round(sd(lit_vec), 2)))

## Comparing Literatuere, MLE, and Bayesian models
print('Comparing Literature, MLE, and Bayesian Models')
model_structure_selection = data.frame()
nll_names = best_model
lit_names = colnames(mod_eval_results)[8:18]
bayes_names = colnames(mod_eval_results)[19:22]
for(i in 1:length(mod_eval_results[ ,1])){
  score_int = min(nll_names[i], na.rm = TRUE)
  best_int = min(mod_eval_results[i,colnames(mod_eval_results) == nll_names], na.rm = TRUE)
  score_lit = min(mod_eval_results[i,8:18], na.rm = TRUE)
  best_lit = lit_names[which(mod_eval_results[i,8:18] == score_lit)]
  score_bayes = min(mod_eval_results[i,19:22], na.rm = TRUE)
  best_bayes = bayes_names[which(mod_eval_results[i,19:22] == score_bayes)]
  best_overall = c('MLE', 'Lit', 'Bayes')[which.min(c(score_int, score_lit, score_bayes))]
  best_mod = c(best_int, best_lit, best_bayes)[which.min(c(score_int, score_lit, score_bayes))]
  write_line = data.frame('best_ll_mod' = best_int, 'score_ensemble' = score_int, 'best_lit_mod' = best_lit, 'score_lit' = score_lit, 'best_bayes_mod' = best_bayes, 'score_bayes' = score_bayes, 'best_model' = best_overall, 'best_mod' = best_mod)
  model_structure_selection = rbind(model_structure_selection, write_line)
}
lit_eval_results_table = aggregate(model_structure_selection$best_lit_mod, by = list(model_structure_selection$best_lit_mod), FUN = length)
best_lit_mod = lit_eval_results_table$Group.1[which.max(lit_eval_results_table$x)]

### Getting summary stats on the best performing literature model
print('Summary Stats of best performing lit mod')
best_lit = mod_eval_results[ ,as.character(best_lit_mod)]
print(paste('range:', round(range(best_lit, na.rm = TRUE), 2)))
best_lit_vec = as.vector(best_lit)
best_lit_vec = best_lit_vec[!is.na(best_lit_vec)]
print(paste('mean:', round(mean(best_lit_vec), 2)))
print(paste('standard deviation:', round(sd(best_lit_vec),2)))

## Getting summary stats for the second best performing literature model
print('Summary Stats of second best performing literature model')
second_best_lit_mod = as.character(lit_eval_results_table$Group.1[order(lit_eval_results_table$x, decreasing = TRUE)[2]])
second_best_lit = mod_eval_results[ ,as.character(second_best_lit_mod)]
second_best_lit_vec = as.vector(second_best_lit)
print(paste('range:', round(range(second_best_lit_vec, na.rm = TRUE), 2)))
print(paste('mean:', round(mean(second_best_lit_vec, na.rm = TRUE), 2)))
print(paste('standard deviation:', round(sd(second_best_lit_vec, na.rm = TRUE), 2)))


## Write results out
save.image(file = file.path(run_results_dir, 'workspace_image_preboot.RData'))

##### Bootstrapping tagging only and prefered models #####
print('Bootstrapping model 5 and prefered model')
boot_iterations = 10000
bootstrap_results = list()

### We'll begin by bootstrapping Model 5 (just tagging data)
## Specifying starting parameters, as well as upper and lower bounds for parameter estimation
#        mu.L, sig.L,  k,  mu.A, sig.A, sig.sci, sig.f, a0, sig.oto, sig.lf
p0 <- c(  70,     1, .10,   1.0,   .10,      1,     0,   0,    0,     0)
lb <- c(  50,   0.1, .05,   0.1,   .05,    0.1,     0,   0,    0,     0)
ub <- c( 110,  15.0, .50,   1.5,   .50,   15.0,     0,   0,    0,     0)

print('Booting Model 5')
# send_push(user = 'uGEHvA4hr37tsrCCtpSv4sUUxVuTqN', message = "model 5")
timer5full = proc.time()
bootstrap_results$booted_param_ests_model5 = bootstrap_growth_params(filename = 'bootstrapped_parameter_estimates_model_5', boot_iterations = boot_iterations, wt.oto = 0, wt.lf = 0, wt.tag = 1, tagdat = tagdat)
bootstrap_results$booted_param_ests_model5withPIFG = bootstrap_growth_params(filename = 'bootstrapped_parameter_estimates_model_5', boot_iterations = boot_iterations, wt.oto = 0, wt.lf = 0, wt.tag = 1, tagdat = tagdat, tagdat2 = tagdat2, wt.tag2 = 1)
bootstrap_results$booted_param_ests_model5justPIFG = bootstrap_growth_params(filename = 'bootstrapped_parameter_estimates_model_5', boot_iterations = boot_iterations, wt.oto = 0, wt.lf = 0, wt.tag = 0, tagdat = tagdat, tagdat2 = tagdat2, wt.tag2 = 1)

boot_time =  (proc.time() -  timer5full)[3] / 60 / 60
# send_push(user = 'uGEHvA4hr37tsrCCtpSv4sUUxVuTqN', message = paste(round(boot_time, digits = 2), "Hours later, bootstrapping model 5  complete!"))

### Now we'll bootstrap the prefered model structure  
## Setting intial params for all data
#        mu.L, sig.L,  k,  mu.A, sig.A, sig.sci, sig.f, a0, sig.oto, sig.lf
p0 <- c(  70,     1, .10,     1,    .1,      1,     0,   0,    1,      1)
lb <- c(  50,   0.1, .05,   0.1,   .05,    0.1,     0, -10,  0.1,    0.1)
ub <- c( 110,  15.0, .50,   1.5,   .50,   15.0,     0,  10,   15,     15)

if (best_integrative_model == 'model 6') {
## 6. Model including all Data sources - Equal weighting to each data type
print('Booting Model 6')
# send_push(user = 'uGEHvA4hr37tsrCCtpSv4sUUxVuTqN', message = "model 6")
timer6 = proc.time()
bootstrap_results$booted_param_ests_model6 = bootstrap_growth_params(filename = 'bootstrapped_parameter_estimates_model_6',boot_iterations = boot_iterations, wt.oto = 1/length(otodat$age), wt.lf = 1/length(lfdat$curr_month_year), wt.tag = 1/dim(tagdat)[1],   otodat = otodat, tagdat = tagdat, pseudolf = pseudo_data)
boot_time =  (proc.time() -  timer6)[3] / 60 / 60
 send_push(user = 'uGEHvA4hr36tsrCCtpSv4sUUxVuTqN', message = paste(round(boot_time, digits = 2), "Hours later, bootstrapping model 6 complete!"))

} else if (best_integrative_model == 'model 7') {
## 7. Model including all Data sources - weighting based on number of sample size
print('Booting Model 7')
# send_push(user = 'uGEHvA4hr37tsrCCtpSv4sUUxVuTqN', message = "model 7")
timer7 = proc.time()
bootstrap_results$booted_param_ests_model7 = bootstrap_growth_params(filename = 'bootstrapped_parameter_estimates_model_7_all_data', boot_iterations = boot_iterations,tagdat=tagdat, otodat=otodat, pseudolf=pseudo_data, wt.oto=1, wt.tag=1, wt.lf=1)
boot_time =  (proc.time() -  timer7)[3] / 60 / 60
 send_push(user = 'uGEHvA4hr37tsrCCtpSv4sUUxVuTqN', message = paste(round(boot_time, digits = 2), "Hours later, bootstrapping model 7 complete!"))

} else if (best_integrative_model == 'model 8') {
## 8. Model including all Data sources treated individually - with equal weighting
print('Booting Model 8')
# send_push(user = 'uGEHvA4hr37tsrCCtpSv4sUUxVuTqN', message = "model 8")
timer8 = proc.time()
bootstrap_results$booted_param_ests_model8 = bootstrap_growth_params(filename = 'bootstrapped_parameter_estimates_model_8_all_data', boot_iterations = boot_iterations, tagdat=tagdat,  otodat=otodat[otodat$source == 'demartini', ], otodat2=otodat[otodat$source == 'ralston and miyamoto', ], otodat3=otodat[otodat$source == 'andrews bomb carbon', ], otodat4=otodat[otodat$source == 'andrews lead radium', ], pseudolf=pseudo_data, pseudolf2=NULL, wt.oto= 1/dim(otodat[otodat$source == 'demartini', ])[1], wt.oto2= 1/dim(otodat[otodat$source == 'ralston and miyamoto', ])[1], wt.oto3=1/dim(otodat[otodat$source == 'andrews bomb carbon', ])[1], wt.oto4=1/dim(otodat[otodat$source == 'andrews lead radium', ])[1], wt.tag = 1/dim(tagdat)[1],  wt.lf = 1/length(pseudolf$curr_month_year), wt.lf2 = 0)
boot_time =  (proc.time() -  timer8)[3] / 60 / 60
 send_push(user = 'uGEHvA4hr37tsrCCtpSv4sUUxVuTqN', message = paste(round(boot_time, digits = 2), "Hours later, bootstrapping model 8 complete!"))

} else if (best_integrative_model== 'model 9') {
## 9. Model including all Data sources treated individually - weighting based on number of sample size
print('Booting Model 9')
# send_push(user = 'uGEHvA4hr37tsrCCtpSv4sUUxVuTqN', message = "model 9")
timer9 = proc.time()
bootstrap_results$booted_param_ests_model9 = bootstrap_growth_params(filename = 'bootstrapped_parameter_estimates_model_9', boot_iterations = boot_iterations, tagdat=tagdat,  otodat=otodat[otodat$source == 'demartini', ], otodat2=otodat[otodat$source == 'ralston and miyamoto', ], otodat3=otodat[otodat$source == 'andrews bomb carbon', ], otodat4=otodat[otodat$source == 'andrews lead radium', ], pseudolf=pseudo_data, pseudolf2=NULL, wt.oto= 1, wt.oto2= 1, wt.oto3=1, wt.oto4=1, wt.tag = 1,  wt.lf = 1, wt.lf2 = 0)
boot_time =  (proc.time() -  timer9)[3] / 60 / 60
 send_push(user = 'uGEHvA4hr37tsrCCtpSv4sUUxVuTqN', message = paste(round(boot_time, digits = 2), "Hours later, bootstrapping model 9 complete!"))

} else if (best_integrative_model == 'model 10') {
## 10. Model without Ralston & Miyamoto - Equal weighting (Because Brett said this was shit!)
print('Booting Model 10')
# send_push(user = 'uGEHvA4hr37tsrCCtpSv4sUUxVuTqN', message = "model 10")
timer10 = proc.time()
bootstrap_results$booted_param_ests_model10 = bootstrap_growth_params(filename = 'bootstrapped_parameter_estimates_model_10', boot_iterations = boot_iterations, tagdat=tagdat, otodat=otodat[otodat$source == 'demartini', ], otodat2=otodat[otodat$source == 'ralston and miyamoto', ], otodat3=otodat[otodat$source == 'andrews bomb carbon', ], otodat4=otodat[otodat$source == 'andrews lead radium', ], pseudolf=pseudo_data, pseudolf2 = NULL, wt.oto= 1/dim(otodat[otodat$source == 'demartini', ])[1], wt.oto2= 0, wt.oto3=1/dim(otodat[otodat$source == 'andrews bomb carbon', ])[1], wt.oto4=1/dim(otodat[otodat$source == 'andrews lead radium', ])[1], wt.tag = 1/dim(tagdat)[1], wt.lf = 1/length(pseudolf$curr_month_year), wt.lf2 = 0)
boot_time =  (proc.time() -  timer10)[3] / 60 / 60
 send_push(user = 'uGEHvA4hr37tsrCCtpSv4sUUxVuTqN', message = paste(round(boot_time, digits = 2), "Hours later, bootstrapping model 10 complete!"))

} else if (best_integrative_model == 'model 11') {
### 11. Model without Ralston & Miyamoto - weighted by n (Because Brett said this was shit!)
print('Booting Model 11')
# send_push(user = 'uGEHvA4hr37tsrCCtpSv4sUUxVuTqN', message = "model 11")
timer11 = proc.time()
bootstrap_results$booted_param_ests_model12 = bootstrap_growth_params(filename = 'bootstrapped_parameter_estimates_model_11', boot_iterations = boot_iterations, tagdat=tagdat, tagdat2 = tagdat2, otodat=otodat[otodat$source == 'demartini', ], otodat2=otodat[otodat$source == 'ralston and miyamoto', ], otodat3=otodat[otodat$source == 'andrews bomb carbon', ], otodat4=otodat[otodat$source == 'andrews lead radium', ], pseudolf=pseudo_data, pseudolf2=NULL, wt.oto= 1, wt.oto2= 0, wt.oto3=1, wt.oto4=1, wt.tag = 1, wt.tag2 = 1, wt.lf = 1, wt.lf2 = 0)
boot_time =  (proc.time() -  timer11)[3] / 60 / 60
 send_push(user = 'uGEHvA4hr37tsrCCtpSv4sUUxVuTqN', message = paste(round(boot_time, digits = 2), "Hours later, bootstrapping model 11 complete!"))
}

## Cleanup
Saving our R environment as an image for easy future reference

save.image(file.path(results_dir, 'Bayesian Workspace.RData'))